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Abstract 

We introduce the quadratic balanced optimization problem (QBOP) which can be 
used to model equitable distribution of resources with pairwise interaction. QBOP is 
strongly NP-hard even if the family of feasible solutions have a very simple structure. 
Several general purpose exact and heuristic algorithms are presented. Results of ex- 
tensive computational experiments are reported using randomly generated quadratic 
knapsack problems as the test bed. These results illustrate the efficacy of our exact and 
heuristic algorithms. We also show that when the cost matrix is specially structured, 
QBOP can be solved as a sequence of linear balanced optimization problems. As a 
consequence, we have several polynomially solvable cases of QBOP. 

Keywords: combinatorial optimization, balanced optimization, knapsack problem, 
bottleneck problems, heuristics. 

1 Introduction 

Let E = {1,2,..., m} be a finite set and J-" be a family of non-empty subsets of E. It is 
assumed that J-' is represented in a compact form of size polynomial in m without explicitly 
listing its elements. For each E E x E, a cost Cij is prescribed. The elements of J-" are 
called feasible solutions and the m x m matrix C = (cjj) is called the cost matrix. Then the 
quadratic balanced optimization problem (QBOP) is to find S E such that 

maxjcjj : {i,j) E S x 8} — minjcjj : {i,j) E S x 3} 

is as small as possible. 
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QBOP is closely related to the balanced optimization problem introduced by Martello et 
al [2T]. To emphasize the difference between the balanced optimization problem of [21J and 
QBOP, we call the former a linear balanced optimization problem (LBOP). Special cases of 
LBOP were studied by many authors [lOl [131 El [13 UHl 1211 [23] . Optimization problems with 
objective functions similar to that of LBOP have been studied by Zeitln [31] for resource 
allocation, by Gupta and Sen [16], Liao and Huang [20], and Tegez and Vlach [291 [301 [31] for 
machine scheduling, by Ahuja [2] for linear programming, and by Scutella [2S] for network 
flows. A generalization of LBOP where elements of E are categorized has been studied by 
Berezny and Lacko [3l [1] and Grineova, Kravecova, and Kulae [15] . Punnen and Nair [26] 
studied LBOP with an additional linear constraint. Punnen and Aneja [25] and Turner et 
al. [32] studied the lexicographic version of LBOP. To the best of our knowledge, QBOP 
has not been studied in literature so far. 

Most of the applications of LBOP discussed in literature translate into applications of 
QBOP by interpreting Cij as the pairwise interaction weight of elements i and j in E. To 
illustrate this, let us consider the following variation of the travel agency example of Martello 
et al. [21]. A North American travel agency is planning to prepare a European tour package. 
Its clients travel from New York to London by a chartered flight. The clients have the option 
to choose a maximum of two tourist locations from an available set 5* of locations. If a client 
chooses locations i and j, then Cij is the total tour time. There are n potential locations 
and the company wishes to choose k = \S\ locations to be included in the package so that 
the duration of tours for any pair of locations is approximately the same. This way, one can 
avoid waiting time of clients in London, after their tour and the whole group can return by 
the same chartered flight. The objective of the tour company can be represented as Mini- 
mizing maxjcjj : (i, j) E S x S} — min{cjj : {i,j) E S x S} while satisfying appropriate 
constraints. Another application of the model arises in balanced portfolio selection for man- 
aging investment accounts where risk estimates on pair of investment opportunities are to 
be considered because of hedging positions. Selection of gift items that are to be packaged 
into pairs of items is yet another application of the model. 

In this paper we study QBOP and propose several general purpose algorithms. The poly- 
nomial solvability of these algorithms are closely related to that of an associated feasibility 
problem. QBOP is observed to be NP-hard even if the family J-" of feasible solutions have 
very simple structure. We also investigate QBOP when the cost matrix C has a decompos- 
able structure, i.e., Cij = + bj or Cij = aibj. In each of these special cases, we show that 
QBOP can be solved in polynomial time whenever the corresponding LBOP can be solved in 
polynomial time. As a consequence, we have O(m^logn) and 0{n^) algorithms for QBOP 
when J-" is chosen as spanning trees of a graph on n nodes and m edges or perfect matchings 
on a Kn,n, respectively. Our general purpose exact algorithms can be modified into heuristic 
algorithms. Some sufficient conditions are derived to speed up these algorithms and their 
effect is analyzed using extensive experimentation in the context of quadratic balanced knap- 
sack problems. We also compared the heuristic solutions with exact solutions and the results 
establish the efficiency of our heuristic algorithms. 

The paper is organized as follows. In section 2 we discuss the complexity of the problem 
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and introduce notations and definitions. Section 3 deals with exact and heuristic algorithms. 
In section 4 we present our we present polynomially solvable special cases. In section 5 we 
discuss the special case of the quadratic balanced knapsack problem. Experimental results 
are presented in section 6 followed by concluding remarks in section 7. 

2 Complexity and notations 

Without loss of generality, we assume that > for otherwise we can add a large constant 
to all Cij values to get an equivalent problem with non-negative cost values. It may be noted 
that when Cij > for all & E x E and Ca = for i & E, QBOP reduces to the 

quadratic bottleneck problem (QBP) [71 |2Z|- QBP is NP-hard even if J-" is the collection of 
all subsets of E with cardinality no more than k for a given k, which depends on m [27j. 
In fact, for such a problem, computing an e-optimal solution is also NP-hard for any e > 
even if Cij G {0, 1} and Ca = [27] . As an immediate consequence, it can be verified that for 
the corresponding instance of QBOP, computing an e-optimal solution is NP-hard for any 
e > 0. In contrast, the corresponding LBOP is polynomially solvable. Thus, the complex- 
ity of QBOP and LBOP are very different and QBOP apparently is a more difficult problem. 

For a given cost matrix C and 5 G J-", we denote 



For a given family of feasible solutions, we use the notation QBOP(C) to indicate that the 
cost matrix under consideration for QBOP is C. Thus, QBOP(C) and QBOP(C*), where 
C ^ C*, are two instances of QBOP with the same family of feasible solutions but different 
cost matrices C and C* respectively. 

For any two real numbers a and (3 such that a < (3 and cost matrix C, let E{C, a, f3) = 



{S* G J-" : Zinin(C, S) > a and Z^g^^{C, S) < /3} and E{C, a, (3) = ■ Cij < a or Cij > f3}. 



Then the quadratic feasibility problem can be stated as follows: "Given two real numbers 
a and (3, where a < /3, test if F{C,a,(3) ^ and produce an 5 G F(C, a,/?) whenever 



Any solution S ^ T can be represented by its incidence vector x = (xi,X2, ■ ■ ■ ,Xm), 
where 



The solution represented by an incidence vector x is denoted by S{x). Let J-huu be the 
convex hull of the incidence vectors of all the elements of J-". Consider the cost matrix C 



^max(C, S) = max{Qj : (i, j) e S x S}, 
^min(C, S) = mm{cij : {i,j) e S X S} and 

Z (C, S) = Zyaa.x{C, S) — Zjam{C, S) . 



F{C,a,(3)^il}. 




1 if 2 G 5, 
otherwise. 



3 



given by 




otherwise. 



1 if {i,j)eE{C, a, (3) 



Then the quadratic feasibihty problem has a 'yes' answer if and only if the optimal objective 
function value of the quadratic programming problem (QPP) 



is zero and if is the corresponding optimal basic feasible solution then S{x^) G F{C, a, 
Thus, the quadratic feasibility problem can be solved by solving the QPP. 

The quadratic feasibility problem can also be viewed as the feasibility version of the linear 
combinatorial optimization problem with conflict pairs (LCOP), where the associated set of 
conflict pairs is precisely E{C, a, i.e, the quadratic feasibility problem has a "yes" answer 
if and only if the set 



is non-empty. For details on the LCOP we refer to PTIISS] . The quadratic feasibility problem 
discussed above is closely related to the quadratic feasibihty problem studied by Punnen and 
Zhang [27] in the context of quadratic bottleneck problems. 

3 Exact and heuristic algorithms for QBOP 

Let us now consider some general results which are used in the subsequent sections to design 
algorithms for QBOP. Let J-'* = {Si, S2, ■ ■ ■ , Sr} be a subset of J-" satisfying the following 
properties: 

(PI) There exists an Si G J-"* which is an optimal solution to QBOP, 

(P2) Zinax(C*, >S'i) < Zinax(C', >S'2) < ■ ■ ■ < •^max ( C, Sr ) . 

For any index k, 1 < k < r, let Tr{k) be the index such that Z^a_^{C, S'^(fc)) — Zmin(C, S'^(fc)) = 
mm{Z^i^^{C, Si) — Zmin(C, Sj) : 1 < i < k}. Then, clearly, STr{r) is an optimal solution to 
QBOP. Let ^7 be a real number such that Q > Z^^in{C, Si) for any optimal solution Si for 
QBOP in T*. 

Theorem 1. For any 1 < k <r, if Z^^^{C, S^^k)) - ^min(C, S^^k)) + ^ < Z^JyC, Sk) then 
Snik) is an optimal solution to QBOP. 



m 



m 




i=l j=l 

subject to X G J-huU 



{x : X E -Fhuib Xi + Xj < 1 for {i,j) G E{C, a, (3)} 



(1) 
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Proof. Suppose S'7r(fc) is not an optimal solution to QBOP. Then there is an optimal solution 
S, to QBOP in J"* such that i > k. Thus, Z^^^{C,Si) - Z„,in{C,Si) < ^max(C, - 
{C,S^(^k))- Then 

(C, Si) < (C, S^(k)) 

< Q + Z^i,x{C, ST^(^k)) ~ Zmm{C, 5'^(fc)) 

^ ■Zmax(C) Sk) ■ 

Thus, by (P2), i < k, a contradiction. □ 

Theorem 2. If S-j^^k) is not an optimal solution to QBOP then there exists an optimal solu- 
tion Sg e T* such that q > k and Z^^^{C, Sg) > Z^J^C, Sk) - ^max(C, + ^mm(C, 5'^(fc)) 

Proof. Since ST,[k) is not optimal, by property (PI) there exists an optimal solution Sg in 
F* such that q> k. Suppose Z^,^{C,S^ < Z^^^{C,Sk) - Z^^^{C, S^(^k)) + Zmm(C, 
Then Z^i^(C, Sg) < Z^^^^^C, Sg) — ^max(C) 'S'7r(fc)) + ■Zmm(C) 'S'7r(fe))- Thus, Zya^^(C, Sg) — 
Zmin{C, Sg) > ^max(C'5 '^^^(fc)) — Z^i^{C, S'^(fe)) establishing that S'^(fc) is also an optimal solu- 
tion to QBOP, a contradiction. □ 

Theorems [1] and [2] assist us in improving the average performance of our algorithms. 
Corresponding results can be obtained by considering solutions that satisfy another set 
of properties. Suppose J-"" = {Si, S2, ■ ■ ■ , S^} be a subset of J-" satisfying the following 
properties. 

(P3) There exists an Si G J-"" which is an optimal solution to QBOP 

(P4) Zinin(C', 'S'l) > Zmiii(C, 5*2) > ■ ■ ■ > Z^i^^C, Sh)- 

Choose an index a{k) such that Zmax(C, Sa(k))-Zmin{C, S„(^k)) = min{Zmax(C, Si)-Zmin{C, Si) : 
1 < i < k}. Then ^^-(/x) is an optimal solution to QBOP. Let A be a real number such that 
A < Z„ia.^{C, Si) for any optimal solution Si for QBOP in 

Theorem 3. For any 1 < k < h, if A - Zrai,^{C, Sa{k)) + ^min(C, S^ik)) > ^mm(C, Sk) then 
Scr{k) is an optimal solution to QBOP. 

Theorem 4. // S'o-(fc) is not an optimal solution to QBOP then there exists an optimal 
solution Sd e 7"° such that d > k and Z^^^{C,Sd) < Z^^^{C, S^(k)) - ^min(C, 5'^(fe)) + 

Zmin{C, Sk) 

The proofs of Theorems |3] and H] can be obtained along the same lines as that of Theo- 
rems [1] and [2] and hence omitted. 

Conditions similar to Theorems [1] to H] have been used by many authors in the context 
of different optimization problems involving linear terms [H [221 IISl 1211 EH]- The effect of 
such conditions are not tested in the context of quadratic type problems. One of the goals 
of our experimental analysis was to test the efficacy of Theorems [1] to H] in the development 
of practical algorithms. 
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3.1 The double threshold algorithm 

The basic idea of this algorithm is similar to that used by Martello et al. |i21j for solving 
LBOP. The difference, however, is that we are using the quadratic feasibility problem 
discussed in the last section instead of a simple feasibility test considered in [21] for LBOP. 
This is a significant deviation as it alters the problem complexity substantially. The validity 
proof of our algorithm follows along the same line as that of the double threshold algorithm 
for LBOP discussed in [21] . We also use the conditions provided in Theorems [1] and [2] to 
enhance our search for an optimal solution. 

Let Wi < W2 < ■ ■ ■ < be an ascending arrangement of distinct elements of the 
cost matrix C and Wp+i = oo. These Wi values are the candidates for .^max(C, S) and 
Zmm{C, S) for any feasible solution 5*. The algorithm performs a bottom-up sequential search 
by maintaining a lower threshold L and an upper threshold U, and tests if F{C, L, U) ^ 
0. If the answer is 'yes', the lower threshold is increased and if the answer is 'no', the 
upper threshold is increased. The lower and upper threshold values are chosen amongst 
{wi,W2, ■ . ■ ,Wp}. At any stage of the algorithm, if a feasible solution is obtained with the 
QBOP objective function value as zero, the algorithm is terminated since we have an optimal 
solution. Let L = W£ and U = Wu for some i and u, £ < u. If F{C,we,Wu) = then U 
is increased to Wu+i- Otherwise, we choose an S* G F{C,W£,Wu) and L can be increased to 
Wv^i, where = minjcjj : (i, j) & S x S}, and the best solution identified so far is updated, 
if necessary. Note that > we. 

We also try to exploit the conditions of Theorems [1] and [2] for early detection of an optimal 
solution or rapid increase in the lower threshold (and hence, possibly the upper threshold). 
Let t < p he the total number of times L is updated and let J-"* = {Si, 5*2, ... , St} be the 
set of solutions generated. The indexes are selected such that Si is generated before Si+i. 
Then J-'* satisfies the properties (PI) and (P2). Thus, the sufficient condition of Theorem [1] 
can be used to detect optimality in any iteration, whenever the condition is satisfied. If 
it is satisfied then the best solution identified so far is indeed optimal and the algorithm 
terminates. Otherwise, we try to increase the lower threshold L (and hence possibly the 
upper threshold U) rapidly using the conditions of Theorem [2l If the algorithm is not 
terminated using any of the conditions discussed above, then the search completes when 
U or L becomes oo and the best solution produced during the search is selected as the 
output which is an optimal solution to QBOP. A formal description of the bottom-up double 
threshold algorithm (BDT algorithm) is given in Algorithm [TJ 

With the assumption that the dominating complexity of this algorithm in each iteration 
is the complexity of testing the condition if F{C,wi,Wu) 7^ or not. If this test can be 
performed in O(0i(m)) time, then the BDT-algorithm terminates in O(m0i(m)) time. For 
most problems of practical interest, testing if F{C,wi,Wu) 7^ or not is NP-hard. By 
performing this test using heuristic algorithms, we can get a heuristic algorithm to solve 
QBOP. The computational issues associated with this approach are discussed in detail in 
the section of experimental analysis of algorithms. 

Just like the BDT algorithm, it is possible to obtain another double threshold algorithm 
using a top-down search. In this case, we start with the upper and lower threshold values at 
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Algorithm 1: The BDT Algorithm 

1 Let wi < W2 < ■ ■ ■ < Wp he an ascending arrangement of distinct values of Cij for 

2 Z ^ 1; u ^ 1; so/ 0; obj -(r- oo; 

3 Compute the parameter ^2; 

4 while I < p and u < p do 



5 a F{C,wi,Wu) then 

6 Choose an S* G F{C,wi,Wu)', 

7 if Z{C, S) < oh] then oh] ^ Z{C, S); sol ^ S; 

8 if ohj = or oh] + VL <wt then return sol; 

9 Choose smallest k such that Wk > max{Zmin(C, S), Z^^i^^{C, S) — ohj}; 

10 I k] 

11 if wi > Wu then u i— k; 

12 else 

13 U ^ U + 1 

14 end if 



15 end while 

16 return ohj and sol 



Wp and systematically decrease the threshold values. The algorithm makes use of Theorems [3] 
and m to improve the search process, along similar lines as in the BDT algorithm where the 
threshold values are systematically increased. The resulting algorithm is called top-down 
double threshold algorithm (TDT algorithm). The detailed description of various steps 
of this algorithm can be easily constructed in view of the BDT algorithm and, therefore, 
omitted. 

3.2 Iterative bottleneck algorithms 

Let us now discuss two additional algorithms for solving QBOP which solve a sequence of 
quadratic bottleneck problems. The worst case complexities of these algorithms , in general, 
are higher than that of the BDT-algorithm, but their average performance is expected to be 
better. A quadratic bottleneck problem of type 1 (QBPl) is defined as 

QBPl: Minimize Z^^^{C,S) 

subject to S E J-'. 

We denote an instance of QBPl with cost matrix C as QBPl(C). The problem QBPl 
was investigated by Burkard [7j and Punnen and Zhang [27j. 

To develop our iterative bottleneck algorithms, we consider a generalization of QBOP, 
where a restriction on the upper threshold is imposed on the feasible solutions. Consider the 
problem 
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QB0P1(C, a): Minimize Z^ax(C, S) - Z^^{C, S) 

subject to S E J^, 

Zmin{C, S) >a 

When a = min{cy : E E x E} QBOPl(C,a) reduces to QBOP(C). Let C be an 
m X m matrix defined by 

^, _ fM if dj < a, 
] Cij otherwise, 

where M is a large number. 

Theorem 5. Let be an optimal solution to QBPl with cost matrix C and q he the index 
such that Wq = Z^i^{C, 5°). 

(1) IfZ^^^{C,S^) = M then QBOPl(C,a) is mfeasible. 

(2) If Z^s^^{C, S^) < M and Zinax(C, 5"^) = Zmin(C, 5*°) then 5° is an optimal solution to 
QBOPlfC, a). 

(3) If conditions (1) and (2) above are not satisfied, then either is an optimal so- 
lution to QBOPl(C,a) or an optimal solution to QBOPl(C','~f) is also optimal to 
QB0P1(C, a), where 7 = Wq+i. 

Proof. The proof of (1) and (2) are straightforward. Let us now prove (3). Let x = {S* G 
: a < Z^i^{C, S) < Zinin(C, S^)}. By definition of x 

Z^UC, 5°) > Z^UC, S) for all S E x- (2) 

Since condition (1) of the theorem is not satisfied, by optimality of to QBPl with cost 
matrix C, we have 

^rnax(C, 5°) < Z^,,{C, S) ioT all S E X- (3) 

Multiply inequality (E]) by — 1 and adding to inequality ([3]) we have Z{C, 5*°) < Z{C, S) for 
all S E X- Thus, either 5° is an optimal solution to QB0P1(C,q;) or there exists an optimal 
solution S to QB0P1(C, a) satisfying Z^i^{C, S) > Z^i^{C, S^) and the result follows. □ 

In view of Theorem [5l we can solve QBOP as a sequence of QBPl problems. In each 
iteration, the algorithm maintains a lower threshold a and constructs a modified cost matrix 
C which depends on the value of a. Then, using an optimal solution to QBPl with cost 
matrix C, the lower threshold is systematically increased until infeasibility with respect 
to the threshold values is reached or optimality of one of the solutions generated so far 
is identified using condition (2) of Theorem [51 Let = {Sp^^, Sp^, . . . , Sp^} be the set of 
solutions generated for various QBPl problems, where Sp^ is generated before S'p.^^, 1 < i < 
t — 1. Then by choosing F* = F^, these solutions satisfy properties (PI) and (P2). Thus, 
Theorem [1] can be used to detect optimality early and Theorem [2] may be used to increase the 
lower threshold rapidly. If the algorithm is not terminated using an optimality condition, 
it compares all the solutions generated by the QBPl solver and outputs the overall best 
solution with respect to the QBOP objective function. The resulting algorithm is called the 
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type 1 iterative bottleneck algorithm (IBl algorithm) and its formal description is given in 
Algorithm [2J 



Algorithm 2: The IBl Algorithm 

1 Let wi < W2 < ■ ■ ■ < Wphe 3ja. ascending arrangement of all distinct values of 

2 C ^ C; obj ^ oo ; so/ ^ 0, M 4— 1 + Wp] Zq = Wp] 

3 Compute the parameter Q; 

4 while Zq ^ M do 

5 Solve QBP1(C"); 

6 if // QBPl(C') is infeasible then return and oo; 

7 else let S be the solution of QBP1(C"); 

8 Zq ^ max{c-j. : (i, j) e S x S}; 

9 if 2;o < ^ then 

10 if Z(C, S) < obj then sol ^ S; obj ^ Z{C, S); 

11 if obj = or obj + Q < Wr-^ then return obj and sol; 

12 L ^ max{Zmm(C, S), Z^JyC, S) - obj}; 

''^^ ^ { M ftherwis^e' '^'^ ^' ^' ^ ^' 

14 end if 

15 end while 

16 return obj and sol 



The IBl algorithm solves at most m problems of the type QBPl. Thus, if QBPl can be 
solved in 0{(j)2{^)) time, then QBOP can be solved in 0{m(j)2{m)) time. By solving QBPl 
using a heuristic, we get a heuristic version of the IBl algorithm. 

QBOP can also be solved as a sequence of quadratic bottleneck problems of the maxmin 
type, which we call a quadratic bottleneck problem of type 2 (QBP2). Formally, QBP2 can 
be stated as follows: 

QBP2: Maximize Z^,^{C,S) 

subject to S ^ T . 

QBP2 can be reformulated as QBPl or the algorithms for QBPl [27J can be modified to 
solve QBP2 directly. 

Now, for any real number /3, consider the problem: 

QB0P2(C, /3): Minimize Z^^^{C, S) - Z^^{C, S) 

subject to S E J^, 

When f3 = max{c^j : e E x E}, QBOP2(C,/3) reduces to QBOP(C). Define the cost 
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matrix C defined by 




—M if Cij > /3 or Cij < a 
Cii otlierwise. 



wliere M is a large number. 

Theorem 6. Let 5*° be an optimal solution to QBP2 with cost matrix C and r he the index 
such that Wr = Z^^xiC, S^). 

1. IfZ^i^{C,S°) = -M then QBOP2(C,/3) is infeasible. 

2. If Z^i^{C, S^) > —M and Z^^^^^C, 5*°) = Z^i^{C, 5°) then is an optimal solution to 
QBOP2(C,(3). 

3. If conditions (1 ) and (2 ) are not satisfied, then either is an optimal solution to 
QB0P2(C,P ) or an optimal solution to QB0P2(C, 'j) is also optimal to QB0P2(C,I3 ) 
where 7 = ■u^+i- 

Tlie proof of tliis tlieorem can be constructed by appropriate modifications in tlie proof 
of Tlieorem O and hence is omitted. 

In view of Theorem [6l we can solve QBOP as a sequence of QBP2 problems. In each iter- 
ation, the algorithm maintains a lower threshold a and an upper threshold /3 and construct a 
modified cost matrix C which depends on the value of a and /3. Using an optimal solution to 
QBP2 with cost matrix C, the upper threshold is systematically decreased and the process 
is continued until infeasibility with respect to the threshold values is reached or optimality 
of one of the solutions generated so far is identified using condition (2) of Theorem [61 Let 

= {Srfi, Sn2, . . . , Sr^^} be the set of solutions generated for various QBP2 problems. The 
indexes are selected such that S*,,- is generated before 5*,;.^^, 1 < i < t — 1. Then by choosing 
^0 _ jp2^ these solutions satisfy properties (P3) and (P4). Thus, Theorem [3] may be used 
to detect optimality early in some cases and Theorem H] may be used to decrease the up- 
per threshold rapidly. If the algorithm is not terminated using an optimality condition, it 
compares the solutions generated by the QBP2 solver and outputs the overall best solution 
with respect to the QBOP objective function. The resulting algorithm is called the type 2 
iterative bottleneck algorithm (IB2-algorithm). A formal description of the IB2 algorithm is 
omitted as it can be obtained by appropriate modifications of the IBl algorithm. 

The IB2 algorithm solves at most m problems of the type QBP2. Thus, if QBP2 can 
also be solved in O(03(m)) time then QBOP can be solved in O(m03(m)) time. By solving 
QBP2 using a heuristic, we get a heuristic version of the IB2 algorithm. 

3.3 The double bottleneck algorithm 

Note that Algorithm IBl sequentially increases the lower threshold value while the algorithm 
IB2 sequentially decreases the upper threshold value. The two algorithms can be combined 
to generate another algorithm that alternately increases the lower threshold and decreases 
the upper threshold. The operations of increasing the lower threshold or decreasing the 
upper threshold are carried out by solving a BQPl and QBP2, respectively. The resulting 
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algorithm is called the double bottleneck algorithm (DB- Algorithm). The validity of the DB- 
algorithm follows from Theorems M and E] and the validity of algorithms IBl and IB2. A 
formal description of the DB Algorithm is given in Algorithm [31 



Algorithm 3: The DB Algorithm 



1 Let wi < W2 < ■ ■ ■ < Wp he an ascending arrangement of all distinct values of 
{dj : EEx E}; 

2 C <— C; oh] ^ oo ; so/ ^ 0, M ^ 1 + Wp] z = Wp] z = Wi, 

3 while z ^ M and z ^ —M do 

4 Solve QBP1(C"); 

5 if // QBPl(C') is infeasible then return and oo; 

6 else let S be the solution of QBF1{C'); 

7 z ^ max{c[j : (z, j) E S x S}; 

8 a zq < M then 

9 if Z{C, S) < obj then sol ^ S; obj ^ Z{C, S); 
10 if obj = then return sol and obj; 



11 4j ^ 



12 end if 



r if > Z^in(C, S), ^^^^ ^ . ^ 

[ M otherwise 



13 Solve QBP2(C"); let S be the resulting solution; 

14 z ^ min{c-j : {i,j) E S x S}; 

15 if 5 > —M then 

16 if Z{C, S) < obj then sol ^ S; obj ^ Z{C, S); 

17 if obj = then return sol and obj; 



18 ^ 



I ""'{^ '['C' for each i,j E E; 

[ —M otherwise 



19 end if 

20 end while 

21 return opt and sol 



4 Polynomially solvable cases 

Let us now consider some special cases of QBOP that can be solved in polynomial time. We 
first consider the decomposable cost matrix, where Cij = + bj for given > 0, 6j > and 
i, j E E. We denote such an instance of a QBOP by QBOP(a + b). Let Z'^{C, S) denote the 
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objective function of QBOP(a + b). Then, 



Z~^{C, S) = maxjcjj : & S x S} — mm{cij : (i, j) & S x S} 

= max{ai + bj : (z, j) E S x S} — minjaj + bj : (i, j) E S x S} 

= maxjoi : i E S} + max{6j : i E S} — minjoi : i E S} — min{6i : i E S} (4) 

= max{aj : i E S} + max{— a, : i E S} + max{6j : i E S} — min{6j : i E S} (5) 

Let Wi be some prescribed weight of i E E and g : T ^ W. Duin and Volgenant fL2\ 
showed 

that combinatorial optimization problems of the type 

COP{g): Minimize max{wj : i E S} + g{S) 

subject to S E J-". 

can be solved in 0{m({m)), where ({m) is the complexity of minimizing g{S) over J-". Note 
that 

Z^{C, S) = maxjaj : i E S} + g{S), 

where 5^(5*) = max{— a-j : i E S} + gi{S) and gi{S) = max{6j : i E S} — min{6j : i E S}. 
But minimizing gi{S) over J-" is precisely the LBOP [21] • Thus, recursively applying the 
results of Duin and Volgenant [121, QBOP(a + b) can be solved in 0{m'^r]{m)) time, where 
ri{m) is the complexity of an LBOP with the same family of feasible solutions as that of the 
QB0P(a + 6). 

The complexity of this algorithm can be improved by careful organization of the compu- 
tations. Let a and {3 be two real numbers such that a < (3. Let Qa,/3 = : a < Oi < 
Let 7 = min{6j : i E Qa,i3} and 6 = max{6j : i E Qa,/?}. Consider the constrained LBOP: 

LBOP (a, P): Minimize max{6j : i E S} — min{6j : i E S} 

subject to s E J^, 

max{6j : i E S} < 6, 
min{6j : i E S} > 'j. 

LB0P(a,/3) can be solved using the double threshold algorithm for LBOP with very sim- 
ple modification by choosing the starting and ending threshold values as 7 and S, respec- 
tively, and the complexity of the algorithm for LBOP in [21j is unaltered when adapted to 
solve LB0P(a,/3). Now we can solve QBOP(a+b) as a sequence of problems of the type 
LB0P(q;,/9). As in [21], we use a double threshold search using the values and simply 
replace the feasibility problem in f2T] by the problem LBOP (a, f3). Thus, if LBOP (a, /3) can 
be solved in 0{ri{m)) time, then QBOP(a + b) can be solved in 0{mr]{m)). 

Note that when feasible solutions are spanning trees of a graph on n nodes and m edges, 
LB0P(a,/3) can be solved in 0(m log n) time [2] and, hence, the resulting QBOP(a + 
b) can be solved in 0(m^ log n) time. When J-" is the family of perfect matchings in a 
complete bipartite graph Kn,n then LBOP (a, /3) can be solved in 0{n^) time [21] and, hence, 
QBOP(a + b) can be solved in O(n^) time. 
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Let us now consider the case when bj for & E X E where Oj > 0, 6j > for 

i = 1,2, ... ,n. The corresponding instance of QBOP is denoted by QB0P(a6). Now, for 
any feasible solution S" G J-", 

Z*{C, S) = max{cjj : E S x S} — min{cy : E S x S} 

= maxjajfej : E S x S} — mm{aibj : E S x S} 

= max{aj : i E S} max{bi : i E S} — min{aj : i E S} min{6j : i E S}. (6) 

By fixing maxjaj : i E S} and minjaj : i E S} between a lower threshold a and an upper 
threshold /3, we can solve QB0P(a6) similar to QBOP (a + 6) as a sequence of problems of 
the type LBOP{a, (3) and the only difference is that we use the objective function Z*{C, S) 
instead of Z^{C,S) when comparing various feasible solutions generated in the course of 
the algorithm. The complexity of the algorithm will be the same as that for the problem 
QBOP(a+b) we discussed earlier. 



5 The quadratic balanced knapsack problem 

Let us now consider a specific case of QBOP called the quadratic balanced knapsack problem 
(QBalKP) which can be defined as follows: 

Minimize maxjcy : EExE,Xi = Xj = l} — min{cij : E E x E,Xi = Xj = 1} 

m 

subject to ttjXj > b (7) 

Xj E {0,1} for J E E. (8) 

By choosing C 2^ such that S E implies J2i£S di > b we can see that QBalKP is an 
instance QBOP. The compact representation of J-" is given by the constraints ([7]) and (|8]). 
QBalKP can be used to model the travel agency example and portfolio selection examples 
discussed in section 1. Since the quadratic bottleneck knapsack problem (QBotKP) [36j is a 
special case of QBalKP and QBotKP is strongly NP-hard, QBalKP is also strongly NP-hard. 
(Note that the LBOP version of QBalKP is solvable in polynomial time.) QBalKP can be 
formulated mixed integer program: 

Minimize u — v 

m 

Subject to ajXj > b, 
i=i 

u > CijVij for E E X E, Cij ^ 0, 

V < CijV.j + M(l - yij) for (i, j) E E x E, dj ^ M, 

Vij ~ Xi <0 for E E X E, 

Uij — Xj <0 for E E X E, 

Xi + Xj - Uij < 1 for (z, j) E E X E, 

Xj E {0, 1}, for J E E, 

0<u,v <M, 
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where M = max{cjj '■ 1 < i, j < n}. Solving the mixed integer programming formulation 
given above becomes difficult as the problem size increases. However, we can use the general 
purpose algorithms developed in the previous section to solve QBalKP. To use the algorithms 
IBl, IB2, and DB, we can make use of the algorithm of Zhang and Punnen [36j as the QBPl 
(QBP2) solver. To apply the BDT and TDT algorithms, we need an algorithm to solve the 
corresponding quadratic feasibility problem. 

Recall that the quadratic feasibility problem is: "Given two real numbers a and where 
a< /3, test if F{C, a, /3) 7^ and produce an 5 e F{C, a, (3) whenever F{C, a, /3) 7^ 0." 

In section 2, we indicated that a quadratic feasibility problem can be solved as a combi- 
natorial optimization problem with conflict pair constraints [HI |35] . We now observe that 
the quadratic feasibility problem for QBKP can be solved by solving the maximum weight 
independent set problem (MWIP) on a graph with node set E and edge set E{a,P). An 
integer programming representation of this MWIP is given below. 



Let optimal solution to the MWIP and z* be its optimal 

objective function value. Then F{C, fi) ^ ^ if and only if z* > b. Thus, we can use 
an MWIP solver to implement the algorithms discussed in the previous section for the 
special case of QBalKP. The solution of the quadratic feasibility problem discussed above is 
closely related to the quadratic feasibility problem studied by Zhang and Punnen [36] for the 
quadratic bottleneck knapsack problem with appropriate differences to handle the QBalKP 
objective. 

6 Computational Experiments 

In this section we report results of extensive experimental analysis conducted on randomly 
generated QBalKP instances. The objective of the experiments is to assess the relative 
performance of various algorithms developed in section 3. We have implemented exact and 
heuristic versions of these algorithms and compared the outcomes in terms of solution quality 
and computational time. Our experiments also examined the effectiveness of the conditions 
provided by Theorems [H [2], [21 and H] for early detection of optimality and rapid advancement 
through the search intervals. All the experiments were conducted on an Intel 17-2600 CPU 
based PC. The algorithms are implemented in C#, and CPLEX 12.4 was used to solve the 
mixed integer programming problems within our implementations. x86-64 instruction set 
was used and concurrency was not allowed in our the algorithms as well as in CPLEX. 

All the algorithms discussed in this paper (except the MIP formulation of QBalKP) 
require a feasibility test procedure and the dominating complexity of these algorithms in 
each iteration is that of this procedure. Recall that a feasibility test answers the question if 
there exists a feasible solution to the QBalKP and, whenever the answer is 'y^s', it generates 



m 



MWIP: 




i=l 

subject to Xi + Xj < 1 for (i, j) G E{a,(3), 
Xj e {0, 1}, for j e E. 
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such a solution. In section 5 we observed that this can be achieved by solving a maximum 
weight independent set problem. We can also use other variations of this approach to test 
feasibility and the empirical behavior of different variations could be different. Since the 
feasibility test is carried out several times in the algorithm, the effect of different variations 
of the feasibility tests could affect the the running time as well as solution quality (for 
heuristic algorithms). For definiteness and simplicity, we have restricted our implementation 
to three different feasibility test procedures which are summarized below: 

FTl: Solve the maximum weight independent set problem as described in section 5 and then 
compare the objective value Xll^i ^i^i "with b. If ^^^^ ajXj > b, the resulting solution 
S{x) is a feasible solution to the QBalKP. Otherwise the answer is 'no'. We used 
CPLEX 12.4 to solve the resulting integer program. 

FT2: Consider the maximum weight independent set problem as described in section 5. 
Choose the objective function coefficients to be zerc0 and add the new constraint 



to the formulation. We used CPLEX 12.4 to solve the resulting constrained maximum 
weight independent set problem. The CPLEX solver stops as soon as a feasible solution 
is found. When using this feasibility test procedure, we set the 'MIP Emphasis' pa- 
rameter of CPLEX to 'Emphasize feasibility over optimality' which in our experiments 
provided the best performance. 

FT3: Solve the integer program defined in FT2 by providing a time limit for the mixed 
integer programming solver. Note that if the solver fails to find a feasible solution in 
the allowed time limit, there is no guarantee that a feasible solution does not exist 
and, thus, using such a procedure in any of the algorithms turns an exact algorithm 
into a heuristic. Because of this heuristic decision, the properties (PI) and (P2) may 
not hold precisely. Nevertheless, we make a heuristic assumption that these properties 
hold and proceed accordingly. 

As indicated earlier, even for these special cases of feasibility tests, the solutions returned 
by these feasibility tests could be different and, thus, not only the execution time of each 
feasibility test but also the optimization process itself for a QBalKP algorithm may vary 
even for the exact feasibility tests FTl and FT2. 

We use the following notations to represent our algorithms under different parameter 
settings. MIP stands for the mixed integer programming formulation of the QBalKP solved 
with CPLEX (the parameters of CPLEX are default). BDT, IB and DB denote the BDT, 
IBl and DB algorithms, respectively, where the effect of Theorem [1] is suppressed. By 
default, we use feasibility test FT2. If '^' is added to the name of an algorithm (such as 

^One can use any objective function in this feasibility test. If the objective function is not a constant, it 
is useful to force the solver to stop after it finds the first feasible solution. However, our experiments showed 
that using the original objective function in this feasibility test slows down the algorithms. 



m 




i=l 
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BDT^ or IB^), feasibility test FTl is used. BDT^ and IB*^ stand for the variations of BDT 
and IB, where early optimality detection is guaranteed by Theorem [T] is enabled. BDT* 
denotes the heuristic version of the BDT algorithm, where t is the time limit prescribed for 
each feasibility test (FT3). IB* and DB* denote the heuristic version of the corresponding 
algorithms, where t is the time limit prescribed for each feasibility test (FT3) within the 
QBPl solving procedure. 

6.1 Test problems 

Since this is the first time when the QBalKP is considered in the literature, we have devel- 
oped a class of test instances for the problem. These test instances are random problems 
constructed as follows. For each test instance, we are given a triple (m, a, s), where m > 1 
is an integer, a > and < s < 1. We first generate an m x m matrix C = (c^j), where 
c'j^j is a normally distributed random integer with mean /i = and standard deviation a as 
given. Then the matrix C = (cij) is generated, where Cij = c'^j — min^sC^s. This guarantees 
that Cij > 0. Then, an m- vector (aj), where at a uniformly distributed random integer in the 
range < aj < 1000, is generated. Finally, b is selected as a uniformly distributed random 
integer in the range [250msJ < b < [750ms \. Observe that E[J2i(zs(^i] = E[b] if S O E 
such that IS"! = ms, where E[x] is the expected value of x. In other words, by varying the 
value of s, one can control the number of non-zeros in an optimal solution to the QBalKP 
instance. 

It may be noted that the instances we generated are symmetric in the following sense; 
replacing Cij with (maxj/j/ Ci'j>)—Cij would not, on average, change the properties of the matrix 
C. Hence, the algorithms IBl and IB2 are expected to show similar average performance. 
Thus, hereafter, we do not discuss the algorithm IB2 and denote IBl as IB. Likewise, BDT 
and TDT algorithms are expected to have similar average performance. Thus, we do not 
consider the TDT algorithm in our experimental analysis and focus on the BDT algorithm. 

Figure [T] indicates relative performance of the BDT and IB algorithms as a function of 
the parameter s. We set m = 100 and a = 100 in this experiment. 

It appears that the instance with either large or small values of s are relatively easy to 
solve. For small values of s, each feasibility test takes only a small amount of time since 
it is easy to find a feasible solution if b is small. For large values of s, each feasibility test 
also takes a relatively small time since many such problems become infeasible. Another 
interesting observation is that the number of iterations of the BDT algorithm almost does 
not depend on s (indeed, even if the problem is infeasible, the BDT algorithm will make p 
iterations) while it varies significantly for the IB algorithm making it significantly faster for 
certain class of instances. 

Figure |2] indicates the performance of the BDT and IB algorithms as a function of the 
parameter a of the instance. The value of s is set in this experiment to s = 0.1 and m = 100. 

One can see that the random instances become harder with increase in the value of a. 
Indeed, a larger a leads to a larger number p of distinct weights Cij which, in turn, increases 
the number of iterations of the algorithm. However, p is limited by and, thus, it grows 
slower than a. In fact, the number of iterations of the BDT algorithm is approximately 
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Figure 1: Analysis of the random instances with different values of parameter s. The hues 
with circle marks indicate the running time of the algorithms and the lines with square marks 
indicate the number of feasibility tests applied within the algorithms. 




100 200 300 400 500 600 700 800 900 1,000 
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Figure 2: Analysis of the random instances with different values of parameter a. The lines 
with circle marks indicate the running time of the algorithms and the lines with square marks 
indicate the number of feasibility tests applied within the algorithms. 
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proportional to p. In contrast, the IB algorithm efficiently handles instances with large p. 
Recall that it needs only O(logp) feasibility tests to solve each QBotKP subproblem 
Hence, for instances with large p, the IB algorithm appears preferable. Also it follows from 
our experiments that the time needed for the feasibility test is almost independent on the 
value of a. 

Based on these preliminary observations, we set s = 0.1 and a = 100 for the rest of the 
experiments to assure that the test problems generated are reasonably hard. 

6.2 Comparison of BDT, IB, DB and MIP 

Let us first evaluate performance of the basic algorithms proposed in this paper. In Table [H 
we report the results of our experiments with the BDT^, BDT, IB^, IB, DB^, DB and MIP 
algorithms. The notations used in various columns of the table are explained below: 

• m is the size of the test instance; 

• 'obj' is the objective value of the optimal solution to the problem; 

• p is the number of distinct Cij values in matrix C; 

• A = maxjj- Cij — minjj- Cij; 

• 'Running time, sec' columns report the running time of each of the algorithms; 

• 'Feasibility tests' columns report the number of times feasibility tests are carried out 
within each of the algorithms. Note that in IB, IB^, DB, and DB^ algorithms, we 
are not explicitly solving feasibility problems. However, feasibility problems of similar 
nature are solved within the QBotKP solver that is used within these algorithms. 
Thus, the number of feasibility tests include the feasibility tests carried out within the 
QBotKP solver used, which is a variation of the algorithm by Zhang and Punnen |36j . 

The last row of the table reports the average values for corresponding columns. However, 
it may be noted that the average running time can not be used to judge the algorithm's 
performance in general because the running times vary significantly from instance to instance. 

The best result (running time and number of tests) for each instance is underlined. 

Feasibility test FT2 provides a better performance than feasibility test FTl in each 
case. Indeed, according to Table [H the optimality of a solution to the maximum weight 
independent set problem does not reduce the number of feasibility tests, while reaching the 
optimal solution clearly takes more time than finding a feasible solution. 

The IB algorithm clearly outperforms all other algorithms for each test instance with 
regards to both the number of feasibility tests and the running time. The MIP algorithm 
turns out to be very slow for any practical instances. 
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Running time, sec Feasibility tests 



m 


obj 


P 


A 




BDT 


IB^ 


IB 


DB^ 


DB 


MIP 


BDT^ 


BDT 


IB^ 


IB 


DB^ 


DB 


50 


64 


453 


689 


2.7 


2.2 


0.6 


(L5 


0.6 


0.5 


1.7 


463 


462 


94 


87 


97 


95 


100 


239 


568 


728 


131.5 


22.0 


56.2 


M 


99.3 


12.2 


7018.1 


582 


581 


128 


125 


235 


209 


150 


204 


622 


772 


476.0 


115.8 


240.1 


37.5 


282.7 


43.1 


13187.2 


642 


643 


191 


188 


227 


228 


200 


181 


672 


778 


938.1 


217.7 


475.3 


38.8 


756.0 


75.1 




684 


684 


117 


101 


207 


174 


250 


206 


701 


807 


3543.5 


980.2 


3806.5 


473.5 


8449.6 


916.4 




719 


720 


182 


178 


347 


348 


Avg. 


179 


603 


755 


1018.4 


267.6 


915.7 


111.7 


1917.6 


209.5 




618 


618 


142 


136 


223 


211 



Table 1: Comparison of basic algorithms. 



Running time, sec Early Iterations Feasibility tests 



m 
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n 


BDT 
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IB 
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IB 
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BDT" 


IB 
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50 


64 


466 
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2.2 
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OA 


/ 


/ 
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432 


10 


9 
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87 


86 
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295 
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8^ 


8.2 


/ 


/ 
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14 


13 
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124 
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346 


115.8 


55.3 


37.5 


34.5 


/ 


/ 


643 


477 


21 


20 


643 


485 


188 


186 
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181 


362 


217.7 


120.1 


38.8 


37.1 


/ 


/ 


684 


497 


12 


11 


684 


505 


101 


99 


250 


206 


348 


980.2 


804.9 


473.5 


503.0 


/ 


/ 


720 


554 


20 


19 


720 


563 


178 


177 


Avg. 


179 


363 


267.6 


200.3 


111.7 


116.6 






618 


485 


15 


14 


618 


494 


136 


134 



Table 2: Analysis of efficiency of the early detection speed up. 



6.3 Early Detection of Optimality 

In this section we report the resuhs of our experiments that assess the effectiveness of The- 
orem [1] for algorithms BDT and IB. In Table El we compare the results of the corresponding 
experiments. The column reports the objective value of the optimal solution to the QBP2 
version of the bottleneck knapsack problem (QBotKP2): 

1} — min{cij : E E x E,Xi = Xj = 1} 

(9) 
(10) 

This is used as the value of the parameter fl in BDT*^ and IB^ algorithms. The 'Early' 
columns report if early detection of optimality happened for the particular algorithm and 
test instance. 'Iterations' columns report the number of iterations of the algorithms. 

In our experiments, early optimality detection happened for every test instance and ev- 
ery algorithm, where such an option was enabled. For the BDT algorithm, early detection 
significantly reduced the number of iterations and improved the overall performance. How- 
ever, it reduces the number of iterations of the IB algorithm only by one in each run. Since 
calculating of fl is approximately equivalent to one iteration of the IB algorithm, the overall 
running time of the procedure almost did not change after enabling early detection. 

Another interesting observation that we can make from Table [2] is that the number of 
iterations of the IB algorithm almost does not depend on the size of the problem. Since each 
such iteration needs O(logp) = O(logm) feasibility tests, the real running time of the IB 
algorithm per iteration is expected to be much higher than that of the BDT algorithm given 
large instances. 

Recall that the value of Q does not need to correspond to an optimal solution of the 
QBotKP2 but may be an upper bound to the optimal objective function value of QBotKP2. 
Observe that the algorithm used in our experiments to solve the QBotKP2 maintains a 
search interval which depends on the values of two indices i, u where i < u. When u — i<l 
the algorithm terminates and guarantees an exact optimal solution [36j. This termination 
condition can be relaxed to u — i < d to obtain a heuristic bound. In our implementation, 
we used the termination condition with u — I < d, where 1 < d < p. Then the resulting 
upper bound Q{d) can be calculated as Q{d) = Wu-i- 

The results of our experiments with the approximate value of Q are reported in Figure [31 
When the value of d is relatively small, the upper bound Q{d) is very close to the exact value 
of Q. For d > 40, the upper bound becomes less accurate and the early detection happens in 
neither DTA^^"') nor IB^'^''^ algorithm. Thus, the effect of using an upper bound n{d) instead 
of the optimal objective value Q is relatively small. 



Maximize minjcjj : G E x E,Xi = Xj 



m 



Subject to > b 



Xj G {0, 1}, for j = 1, 2, . . . , m. 
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Figure 3: Analysis of the early optimality detection for different values of d. Observe that 
BDT^(i) = BDT^, and for = p we replace BDT^(p) and IB^(p) with the BDT and IB 
algorithms, respectively. The relative running time is calculated as td/ti, where is the 
running time of the BBT*^*^"^) or IB'^^'^^ algorithm and ti is the running time of the BDT or 
IB algorithm, respectively. The experiments were conducted for a random instance of size 
m = 150. 

6.4 Comparison of Heuristics 

Recall that, by setting a time limitation to CPLEX when running a feasibility test, one can 
speed up the BDT, IB and DB algorithms at the cost of loosing optimality guarantee. In 
Figure HI we show how the solution quality depends on the running time for each of the 
BDT*, IB* and DB* algorithms. 

In our experiments, the running time and the solution quality of each of the heuristics 
monotonically depend on the parameter t. Thus, the balance between the solution quality 
and the running time in the proposed heuristics can be efficiently controlled by t. Note that 
the IB* algorithm, likewise its exact version, shows the best performance among the proposed 
heuristics. 

We also compared the BDT*, IB* and DB* heuristics with the MIP algorithm given a 
time limitation (in this case, we set the 'MIP Emphasis' parameter of CPLEX to 'Emphasize 
feasibility over optimality'). However, in our experiments, such a MIP heuristic showed very 
poor performance. 

7 Conclusion 

We introduced the combinatorial optimization model QBOP which can be used to model 
equitable distribution problems with pairwise interactions. The problem is strongly NP-hard 
even if the family of feasible solutions have a simple structure such as the collection of all 
subsets of a finite set with an upper bound on the cardinality of these subsets. Several exact 
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Running time, sec 

Figure 4: Performance of the BDT*, IB* and DB* heuristics. The experiment is conducted 
for a random instance of size m — 500. The values of t, in seconds, are reported near each 
node. 

and heuristic algorithms are provided along with detailed experimental analysis in the case 
of quadratic knapsack problems. Special cases of the problem with decomposable type cost 
matrices are discussed. It is shown that the complexity of the resulting QBOP depends on 
that of the corresponding LBOP. 

By exploiting special structure of the family of feasible solutions J-" and the structure of 
C, one may be able to obtain improved algorithms. These are interesting topics to investigate 
further, especially when real life situations warrant the study of such problems. 
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